BTE Appendix Table and Figures
Zeyi Han
5/31/2021
- Data
- 1. Appendix I Table S1. Gamm models for Density Data
- 2. Appendix I Table S2. Quantifying temperature and nutrient effects on species temporal ecological and phenotypic dynamics using ARMA GAMM models.
- 3. Appendix I Table S3. Linear Models for Ecological Dynamics Descriptor
- 4. Appendix II Figure S1-14 CCM convergence plot for each species pair
- 5. Appendix II Table S1 Mean CCM Prediction at Largest Library Size by Temperature
- 6. Appendix III Table S1 All Strong Interaction in the FW
Data
Data read in
We used log density to standardize the variance in protists density time series. Since in some days, the protists density is 0, and will result -INF by log(), we added 1 to all the density before calculating the log density for model fitting. We minused the 1 after exponentiating the model predictions.
### data read in------
## density
milou_dens0 <- read.csv("/Users/zeyi/Documents/phd/Projects/Milou/02_Clean_data/Density.csv")
milou_dens <- milou_dens0 %>% mutate(tre=as.factor(tre), temp = as.factor(temp), nut = as.factor(nut), all_rep = rep(1:24, time = 24), log_density = log(density + 1)) %>% select(-1)
tetra <- milou_dens[which(milou_dens$Class=="Tetrahymena"),]
eup <- milou_dens[which(milou_dens$Class=="Euplotes"),]
bac0 <- read.csv("/Users/zeyi/Documents/phd/Projects/Milou/02_Clean_data/Milou_bac_OD.csv")
bac1 <- bac0 %>% mutate(tre = as.factor(tre), temp = as.factor(temp), nut = as.factor(nut), all_rep = rep(1:24, time = 11), temp_num = as.numeric(rep(c(22,25), each = 6, times = 22)), nut_num = as.numeric(rep(c(1,0.5), each = 12, times = 11)))We calculated the mean and standard variation for bacteria OD after day 5. One data point that is more than 3 standard deviation away from the mean is removed.
sd <- sd(bac1$OD[which(bac1$day > 5)], na.rm = TRUE)
mean <- mean(bac1$OD[which(bac1$day > 5)], na.rm = TRUE)
bac <- bac1[-which(bac1$OD > mean+3*sd),]Since after population collapse, Tetrahymena are functionally extinct and so, even if we occasionally found individual, trait distributions are unreliable. For traits observation made after each Tetrahymena treatment collapsed (once population density became 0 ind/ml) that with less than 10 number of individual, we considered them unreliable and, thus, removed for future analysis.
## trait
t_summ0 <- fread("/Users/zeyi/Documents/phd/Projects/Milou/02_Clean_data/Milou_trait_summary.csv")
t_summ <- t_summ0 %>% mutate(rep = as.character(rep),
tre = as.factor(tre),temp_num = as.numeric(temp_num), nut_num = as.numeric(nut_num),rep = as.factor(rep), tre = as.factor(tre), all_rep = rep(1:24, time = 24), temp = as.factor(temp), nut = as.factor(nut))
tetraSumm <- t_summ %>% filter(Class =="Tetrahymena")
eupSummtrait <- t_summ %>% filter(Class =="Euplotes")
# replace NA for tetra trait}
# Half 22C
tetraSumm[, 14:20][which(tetraSumm$tre=="Half22C" & tetraSumm$day>10 & tetraSumm$num_ind<10)] <- NA
# Full 22C
tetraSumm[, 14:20][which(tetraSumm$tre=="Full22C" & tetraSumm$day>11 & tetraSumm$num_ind<10)] <- NA
# Half 25C
tetraSumm[, 14:20][which(tetraSumm$tre=="Half25C" & tetraSumm$day>8 & tetraSumm$num_ind<10)] <- NA
# Full 25C
tetraSumm[, 14:20][which(tetraSumm$tre=="Full25C" & tetraSumm$day>8 & tetraSumm$num_ind<10)] <- NA
# trait group mean
tetraSummMean <- tetraSumm %>% na.omit %>% group_by(Class,temp,nut,tre,day) %>% summarise(mArea = mean(area)) %>% filter(day <= 9)
eupSummMean <- eupSummtrait %>% na.omit %>% group_by(Class,temp,nut,tre,day) %>% summarise(mArea = mean(area)) %>% filter(day <= 9) %>% filter(day <= 9)
# cropped time series for second version of CCM
tetraSumm_crop <- tetraSumm %>% filter(day <= 9)
eupSumm_crop <- eupSummtrait %>% filter(day <=9)Calculate mean density of Bacteria, Tetrahymena, and Euplotes
bac.2 <- bac %>%
rename(density = OD) %>%
select(tre:density,treatNum) %>%
mutate(temp = as.factor(temp))
teb.all <- milou_dens %>%
select(tre:density,temp_num,-temp,treatNum) %>%
rename(temp = temp_num) %>%
mutate(temp = as.factor(temp), nut = as.factor(nut),day = as.numeric(day))%>%
full_join(bac.2) %>%
na.omit() %>%
group_by(Class,tre,treatNum, nut, temp,day) %>%
summarise(mean.d = mean(density),
log.d = log(mean.d +1)) %>%
mutate(day = as.numeric(day))
teb.mean <- teb.all %>%
select(Class:mean.d) %>%
spread(key=Class, value = mean.d)1. Appendix I Table S1. Gamm models for Density Data
We used GAMM model to select for the detection of temporal autocorrelation and correctly accounting for that. To do so, we treated replicates as random effects and day and treatment as fixed effects. Temporal autocorrelation was accounted for using Autoregressive Moving Average (ARMA) correlation structure (of orders 1, 2 and 3) using “corARMA” in R package “nlme” to apply in the GAMM models using the ‘nlme’ R package (version 3.1-148). GAMM models smoothed by treatment more closely reproduce the observed dynamics and are thus better suited for the detection autocorrelation in the time series data, as well as for time series analysis using CCM.
1.1 Bacteria Density ARMA selection
AR = 2 is the most suitable for Bacteria Density.
#bac1 <- bac1 %>% mutate(tre = as.factor(tre), temp_num = as.numeric(temp), nut_num = as.numeric(rep(c(1,0.5), each = 12, time = 11 )))
bm0 <- gamm(OD ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = bac)
summary(bm0$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## OD ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0251377 0.0009284 27.077 < 2e-16 ***
## treFull25C -0.0010281 0.0013129 -0.783 0.435
## treHalf22C -0.0074084 0.0013129 -5.643 5.71e-08 ***
## treHalf25C -0.0131807 0.0013196 -9.989 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.103 5.103 20.433 < 2e-16 ***
## s(day):treFull25C 3.915 3.915 17.571 3.57e-12 ***
## s(day):treHalf22C 3.775 3.775 7.652 2.55e-05 ***
## s(day):treHalf25C 3.484 3.484 3.269 0.0118 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.611
## Scale est. = 4.6527e-05 n = 219bm1 <- gamm(OD ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 1))
summary(bm1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## OD ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0251262 0.0008449 29.738 < 2e-16 ***
## treFull25C -0.0008859 0.0011946 -0.742 0.459
## treHalf22C -0.0073185 0.0011946 -6.126 4.73e-09 ***
## treHalf25C -0.0130901 0.0012017 -10.893 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 4.944 4.944 23.480 < 2e-16 ***
## s(day):treFull25C 3.807 3.807 20.443 1.48e-13 ***
## s(day):treHalf22C 3.795 3.795 8.333 8.12e-06 ***
## s(day):treHalf25C 3.449 3.449 4.232 0.00991 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.607
## Scale est. = 4.6739e-05 n = 219#AR2
bm2 <- gamm(OD ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 2))
summary(bm2$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## OD ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.025277 0.001201 21.046 < 2e-16 ***
## treFull25C -0.001070 0.001698 -0.630 0.529
## treHalf22C -0.007480 0.001698 -4.405 1.74e-05 ***
## treHalf25C -0.013237 0.001703 -7.774 4.18e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 6.065 6.065 17.671 < 2e-16 ***
## s(day):treFull25C 4.793 4.793 14.396 1.71e-11 ***
## s(day):treHalf22C 4.098 4.098 8.088 4.19e-06 ***
## s(day):treHalf25C 3.636 3.636 4.073 0.00935 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.616
## Scale est. = 4.6931e-05 n = 219#AR3
bm3 <- gamm(OD ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 3))
summary(bm3$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## OD ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.025270 0.001204 20.981 < 2e-16 ***
## treFull25C -0.001055 0.001703 -0.620 0.536
## treHalf22C -0.007468 0.001703 -4.386 1.88e-05 ***
## treHalf25C -0.013224 0.001708 -7.744 5.00e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 6.036 6.036 17.802 < 2e-16 ***
## s(day):treFull25C 4.794 4.794 14.515 1.36e-11 ***
## s(day):treHalf22C 4.112 4.112 8.203 3.41e-06 ***
## s(day):treHalf25C 3.654 3.654 4.108 0.00895 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.616
## Scale est. = 4.6821e-05 n = 219## df AIC
## bm1$lme 15 -1489.672
## bm2$lme 16 -1499.449
## bm3$lme 17 -1497.476
## bm0$lme 14 -1489.6331.2 Tetrahymena Density ARMA selection
We tested gamm models on Tetrahymena as we tested on bacteria. AR = 1 has the lowest AIC.
# no autocorrelation accounted for
tm0 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = tetra)
summary(tm0$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8942 0.1343 36.450 < 2e-16 ***
## treFull25C -1.3444 0.1305 -10.301 < 2e-16 ***
## treHalf22C -0.4596 0.1305 -3.522 0.000508 ***
## treHalf25C -1.6974 0.1305 -13.006 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 7.630 7.630 180.0 <2e-16 ***
## s(day):treFull25C 8.225 8.225 186.5 <2e-16 ***
## s(day):treHalf22C 7.525 7.525 167.7 <2e-16 ***
## s(day):treHalf25C 7.657 7.657 193.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.60466 n = 288#AR1
tm1 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 1))
summary(tm1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8844 0.1479 33.025 < 2e-16 ***
## treFull25C -1.3269 0.1846 -7.190 7.4e-12 ***
## treHalf22C -0.4505 0.1845 -2.442 0.0153 *
## treHalf25C -1.6850 0.1845 -9.132 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 8.063 8.063 116.0 <2e-16 ***
## s(day):treFull25C 8.561 8.561 110.7 <2e-16 ***
## s(day):treHalf22C 7.973 7.973 104.0 <2e-16 ***
## s(day):treHalf25C 8.155 8.155 102.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.59689 n = 288#AR2
tm2 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 2))
summary(tm2$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8780 0.1433 34.050 < 2e-16 ***
## treFull25C -1.3148 0.1734 -7.581 6.61e-13 ***
## treHalf22C -0.4449 0.1734 -2.566 0.0109 *
## treHalf25C -1.6750 0.1734 -9.661 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 8.015 8.015 119.4 <2e-16 ***
## s(day):treFull25C 8.546 8.546 116.6 <2e-16 ***
## s(day):treHalf22C 7.911 7.911 108.6 <2e-16 ***
## s(day):treHalf25C 8.127 8.127 107.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.59936 n = 288#AR3
tm3 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 3))
summary(tm3$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.8784 0.1427 34.187 < 2e-16 ***
## treFull25C -1.3156 0.1724 -7.632 4.81e-13 ***
## treHalf22C -0.4449 0.1723 -2.582 0.0104 *
## treHalf25C -1.6758 0.1723 -9.724 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 8.021 8.021 119.9 <2e-16 ***
## s(day):treFull25C 8.548 8.548 117.0 <2e-16 ***
## s(day):treHalf22C 7.917 7.917 109.1 <2e-16 ***
## s(day):treHalf25C 8.132 8.132 107.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.59929 n = 288## df AIC
## tm0$lme 14 833.3562
## tm1$lme 15 779.9543
## tm2$lme 16 779.4673
## tm3$lme 17 781.44761.3 Euplotes Density ARMA selection
Same models as above. Gamm model without AR structure has the lowest AIC.
eup$log_density <- log(eup$density + 1)
em0 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = eup)
summary(em0$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.10965 0.05723 54.332 < 2e-16 ***
## treFull25C 0.52813 0.08094 6.525 3.56e-10 ***
## treHalf22C 0.18298 0.08094 2.261 0.0246 *
## treHalf25C 0.51999 0.08094 6.424 6.31e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.076 5.076 269.7 <2e-16 ***
## s(day):treFull25C 7.689 7.689 188.7 <2e-16 ***
## s(day):treHalf22C 6.720 6.720 174.7 <2e-16 ***
## s(day):treHalf25C 5.424 5.424 226.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.23258 n = 288#AR1
em1 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = eup, correlation = corARMA(form = ~ day|tre/rep, p = 1))
summary(em1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.10916 0.05902 52.679 < 2e-16 ***
## treFull25C 0.52886 0.08347 6.336 1.04e-09 ***
## treHalf22C 0.18273 0.08347 2.189 0.0295 *
## treHalf25C 0.52155 0.08347 6.248 1.69e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.086 5.086 253.0 <2e-16 ***
## s(day):treFull25C 7.725 7.725 176.7 <2e-16 ***
## s(day):treHalf22C 6.743 6.743 163.8 <2e-16 ***
## s(day):treHalf25C 5.437 5.437 212.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.23285 n = 288#AR2
em2 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = eup, correlation = corARMA(form = ~ day|tre/rep, p = 2))
summary(em2$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.10765 0.06325 49.131 < 2e-16 ***
## treFull25C 0.52931 0.08946 5.916 1.04e-08 ***
## treHalf22C 0.18532 0.08946 2.071 0.0393 *
## treHalf25C 0.52230 0.08945 5.839 1.57e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.181 5.181 221.6 <2e-16 ***
## s(day):treFull25C 7.868 7.868 153.6 <2e-16 ***
## s(day):treHalf22C 6.899 6.899 141.3 <2e-16 ***
## s(day):treHalf25C 5.531 5.531 186.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.23224 n = 288#AR3
em3 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = eup, correlation = corARMA(form = ~ day|tre/rep, p = 3))
summary(em3$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## log_density ~ s(day, by = tre, k = 10) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.10569 0.06136 50.610 < 2e-16 ***
## treFull25C 0.53119 0.08684 6.117 3.51e-09 ***
## treHalf22C 0.18713 0.08682 2.155 0.0321 *
## treHalf25C 0.52492 0.08679 6.048 5.11e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.159 5.159 230.7 <2e-16 ***
## s(day):treFull25C 7.921 7.921 157.5 <2e-16 ***
## s(day):treHalf22C 6.903 6.903 146.1 <2e-16 ***
## s(day):treHalf25C 5.502 5.502 191.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.948
## Scale est. = 0.23251 n = 288## Model df AIC BIC logLik Test L.Ratio p-value
## em0$lme 1 14 512.2403 563.5218 -242.1202
## em1$lme 2 15 513.8043 568.7487 -241.9021 1 vs 2 0.4360305 0.5090
## em2$lme 3 16 514.1925 572.7999 -241.0963 2 vs 3 1.6117600 0.2042
## em3$lme 4 17 515.6682 577.9386 -240.8341 3 vs 4 0.5242958 0.46901.4 Tetrahymena Area ARMA selection
Same models as above but for Tetrahymena phenotypic data. Gamm model with AR = 3 structure has the lowest AIC.
# model with interactive effect from 4 treatments
# No AR
tt0 <- gamm(area ~ s(day, by = tre, k =7) + tre ,random = list(rep= ~ 1), data = tetraSumm)
summary(tt0$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 7) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 789.378 10.063 78.441 < 2e-16 ***
## treFull25C -2.988 25.637 -0.117 0.90737
## treHalf22C -40.946 14.298 -2.864 0.00476 **
## treHalf25C -131.137 31.249 -4.197 4.51e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.244 5.244 46.56 <2e-16 ***
## s(day):treFull25C 3.988 3.988 46.53 <2e-16 ***
## s(day):treHalf22C 5.184 5.184 35.51 <2e-16 ***
## s(day):treHalf25C 3.754 3.754 30.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.808
## Scale est. = 4928.9 n = 180# AR = 1
tt1 <- gamm(area ~ s(day, by = tre, k =7) + tre,random = list(rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|tre/rep, p = 1))
summary(tt1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 7) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 791.929 8.094 97.847 < 2e-16 ***
## treFull25C -4.794 24.169 -0.198 0.843030
## treHalf22C -42.181 11.422 -3.693 0.000305 ***
## treHalf25C -131.298 29.060 -4.518 1.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.075 5.075 66.67 <2e-16 ***
## s(day):treFull25C 3.931 3.931 64.61 <2e-16 ***
## s(day):treHalf22C 4.961 4.961 54.14 <2e-16 ***
## s(day):treHalf25C 3.686 3.686 39.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.803
## Scale est. = 4922.3 n = 180# AR = 2
tt2 <- gamm(area ~ s(day, by = tre, k =7) + tre,random = list(rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|tre/rep, p = 2))
summary(tt2$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 7) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 792.868 8.386 94.544 < 2e-16 ***
## treFull25C -4.129 24.619 -0.168 0.86702
## treHalf22C -42.913 11.952 -3.590 0.00044 ***
## treHalf25C -132.407 29.271 -4.523 1.19e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.186 5.186 61.91 <2e-16 ***
## s(day):treFull25C 3.984 3.984 62.04 <2e-16 ***
## s(day):treHalf22C 5.111 5.111 48.94 <2e-16 ***
## s(day):treHalf25C 3.725 3.725 39.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.804
## Scale est. = 4883.7 n = 180# AR = 3
tt3<- gamm(area ~ s(day, by = tre, k =7) + tre,random = list(rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|tre/rep, p = 3))
summary(tt3$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 7) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 790.886 10.898 72.573 < 2e-16 ***
## treFull25C -9.502 25.069 -0.379 0.70517
## treHalf22C -40.700 15.511 -2.624 0.00954 **
## treHalf25C -128.130 31.031 -4.129 5.87e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 5.038 5.038 76.78 <2e-16 ***
## s(day):treFull25C 3.842 3.842 108.01 <2e-16 ***
## s(day):treHalf22C 4.667 4.667 67.48 <2e-16 ***
## s(day):treHalf25C 3.708 3.708 55.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.796
## Scale est. = 5290.6 n = 180## df AIC
## tt0$lme 14 2139.415
## tt1$lme 15 2131.184
## tt2$lme 16 2132.748
## tt3$lme 17 2126.7371.5 Euplotes Area ARMA Selection
Same models as above but for Euplotes phenotypic data. Gamm model with AR = 1 structure has the lowest AIC.
# interactive effects from different treatment
# No AR
et0 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSummtrait)
summary(et0$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 9) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7885.4 133.6 59.020 < 2e-16 ***
## treFull25C -567.7 186.1 -3.051 0.00253 **
## treHalf22C -237.9 186.6 -1.275 0.20356
## treHalf25C -801.4 186.0 -4.308 2.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 3.991 3.991 18.58 1.12e-13 ***
## s(day):treFull25C 4.501 4.501 14.09 3.55e-11 ***
## s(day):treHalf22C 4.750 4.750 18.50 4.39e-15 ***
## s(day):treHalf25C 4.893 4.893 15.02 6.80e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.543
## Scale est. = 1.1184e+06 n = 267# AR = 1
et1 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|tre/rep, p = 1))
summary(et1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 9) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7868.6 119.8 65.696 < 2e-16 ***
## treFull25C -555.9 165.4 -3.360 0.000903 ***
## treHalf22C -192.1 165.9 -1.158 0.248096
## treHalf25C -779.3 165.6 -4.707 4.21e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 4.080 4.080 21.45 9.03e-16 ***
## s(day):treFull25C 4.587 4.587 16.65 3.04e-13 ***
## s(day):treHalf22C 4.930 4.930 22.57 < 2e-16 ***
## s(day):treHalf25C 4.962 4.962 18.44 5.91e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.544
## Scale est. = 1.1033e+06 n = 267# AR = 2
et2 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|tre/rep, p = 2))
summary(et2$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 9) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7872.2 118.2 66.587 < 2e-16 ***
## treFull25C -559.3 163.1 -3.428 0.000712 ***
## treHalf22C -199.0 163.7 -1.216 0.225286
## treHalf25C -783.9 163.3 -4.801 2.75e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 4.064 4.064 21.72 6.31e-16 ***
## s(day):treFull25C 4.568 4.568 16.66 3.26e-13 ***
## s(day):treHalf22C 4.908 4.908 22.68 < 2e-16 ***
## s(day):treHalf25C 4.953 4.953 18.83 4.06e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.544
## Scale est. = 1.1038e+06 n = 267# AR = 3
et3 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|tre/rep, p = 3))
summary(et3$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = tre, k = 9) + tre
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7874.9 119.8 65.723 < 2e-16 ***
## treFull25C -560.5 165.6 -3.384 0.000831 ***
## treHalf22C -199.6 166.1 -1.202 0.230679
## treHalf25C -784.8 165.7 -4.737 3.69e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):treFull22C 4.065 4.065 21.10 1.7e-15 ***
## s(day):treFull25C 4.593 4.593 16.95 1.7e-13 ***
## s(day):treHalf22C 4.913 4.913 23.17 < 2e-16 ***
## s(day):treHalf25C 4.972 4.972 19.06 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.544
## Scale est. = 1.105e+06 n = 267## df AIC
## et0$lme 14 4559.905
## et1$lme 15 4556.980
## et2$lme 16 4558.873
## et3$lme 17 4560.7672. Appendix I Table S2. Quantifying temperature and nutrient effects on species temporal ecological and phenotypic dynamics using ARMA GAMM models.
Because we have only two temperature and nutrients levels, we treated them as categorical variables. Here we present model comparison using AICc on models that accounts for all combinations of additive and/or interactive effects temperature and nutrient. We apply the ARMA structure selected in these models.
2.1 Bacteria Density
Model 7 with temperature additive and nutrient interactive effects (nutrients within smooth term) is the most parsimonious model.
# 1. Temp additive effect
b.od1 <- gamm(OD ~ s(day, k = 10) + temp, random = list(all_rep = ~1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od1$gam)
# 2. Nut additive effect
b.od2 <- gamm(OD ~ s(day, k = 10) + nut,random = list(all_rep= ~ 1),correlation = corARMA(form = ~ day|all_rep, p = 2), data = bac)
summary(b.od2$gam) # nutrient not significant
# 3. Temp and Nut additive effects
b.od3 <- gamm(OD ~ s(day, k = 10) + nut + temp,random = list(all_rep= ~ 1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od3$gam) # nutrient not significant
# 4. Temp interactive effect ONLY (Temp within smooth term)
b.od4 <- gamm(OD ~ s(day, by = temp, k = 10) + temp, random = list(all_rep = ~1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od4$gam)
# 5. Nut interactive effect ONLY (Nut within smooth term)
b.od5 <- gamm(OD ~ s(day, by = nut, k = 10) + nut , random = list(all_rep = ~1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od5$gam)
# 6. Temp interactive + nut additive (Temp within smooth term)
b.od6 <- gamm(OD ~ s(day, by = temp, k = 10) + nut+ temp, random = list(all_rep = ~1),correlation = corARMA(form = ~ day|all_rep, p = 2), data = bac)
summary(b.od6$gam)
# 7. Nut interactive + temp additive (Nut within smooth term)
b.od7 <- gamm(OD ~ s(day, by = nut, k = 10) + temp+ nut, random = list(all_rep = ~1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od7$gam)
# 8. nut * temp outside of smooth term.
b.od8 <- gamm(OD ~ s(day, k = 10) + temp*nut, random = list(all_rep = ~1), correlation = corARMA(form = ~ day|all_rep, p = 2),data = bac)
summary(b.od8$gam)
# Model Comparison
AIC(b.od1$lme, b.od2$lme,b.od3$lme, b.od4$lme,b.od5$lme,b.od6$lme,b.od7$lme, b.od8$lme)2.2 Tetrahymena Density
Model 6 with nutrient additive and temperature interactive effects (temperature within smooth term) is the most parsimonious model.
# 1. Temp additive effect
t.dens1 <- gamm(log_density ~ s(day, k = 10) + temp, random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|temp/all_rep, p = 1))
summary(t.dens1$gam)
# 2. Nut additive effect
t.dens2 <- gamm(log_density ~ s(day, k = 10) + nut,random = list(all_rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|nut/all_rep, p = 1))
summary(t.dens2$gam)
# 3. Temp and Nut additive effects
t.dens3 <- gamm(log_density ~ s(day, k = 10) + nut + temp,random = list(all_rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(t.dens3$gam)
# 4. Temp interactive effect ONLY (Temp within smooth term)
t.dens4 <- gamm(log_density ~ s(day, by = temp, k = 10) + temp , random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|temp/all_rep, p = 1))
summary(t.dens4$gam)
# 5. Nut interactive effect ONLY (Nut within smooth term)
t.dens5 <- gamm(log_density ~ s(day, by = nut, k = 10) + nut, random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|nut/all_rep, p = 1))
summary(t.dens5$gam)
# 6. Temp interactive + nut additive (Temp within smooth term)
t.dens6 <- gamm(log_density ~ s(day, by = temp, k = 10) + nut + temp, random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(t.dens6$gam)
# 7. Nut interactive + temp additive (Nut within smooth term)
t.dens7 <- gamm(log_density ~ s(day, by = nut, k = 10) + temp + nut, random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(t.dens7$gam)
# 8. temp*nut outside of smooth term
t.dens8 <- gamm(log_density ~ s(day, k = 10) + temp*nut, random = list(all_rep = ~1), data = tetra, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(t.dens8$gam)
# model comparison
AIC(tm1$lme,t.dens1$lme, t.dens2$lme,t.dens3$lme, t.dens4$lme,t.dens5$lme, t.dens6$lme,t.dens7$lme,t.dens8$lme)2.3 Euplotes Density
Model 4 with temperature interactive effects only (temperature within smooth term) is the most parsimonious model.
# 1. Temp additive effect
e.dens1 <- gamm(log_density ~ s(day, k = 10) + temp, random = list(all_rep = ~1), data = eup)
summary(e.dens1$gam)
# 2. Nut additive effect
e.dens2 <- gamm(log_density ~ s(day, k = 10) + nut,random = list(all_rep= ~ 1), data = eup)
summary(e.dens2$gam) # nutrient not significant
# 3. Temp and Nut additive effects
e.dens3 <- gamm(log_density ~ s(day, k = 10) + nut + temp,random = list(all_rep= ~ 1), data = eup)
summary(e.dens3$gam) # nutrient not significant
# 4. Temp interactive effect ONLY (Temp within smooth term)
e.dens4 <- gamm(log_density ~ s(day, by = temp, k = 10) + temp, random = list(all_rep = ~1), data = eup)
summary(e.dens4$gam)
# 5. Nut interactive effect ONLY (Nut within smooth term)
e.dens5 <- gamm(log_density ~ s(day, by = nut, k = 10) + nut, random = list(all_rep = ~1), data = eup)
summary(e.dens5$gam)
# 6. Temp interactive + nut additive (Temp within smooth term)
e.dens6 <- gamm(log_density ~ s(day, by = temp, k = 10) + nut + temp, random = list(all_rep = ~1), data = eup)
summary(e.dens6$gam) # nutrient is not significant, temperature is
# 7. Nut interactive + temp additive (Nut within smooth term)
e.dens7 <- gamm(log_density ~ s(day, by = nut, k = 10) + temp+ nut, random = list(all_rep = ~1), data = eup)
summary(e.dens7$gam)
# 8. Nut * temp outside of smooth
e.dens8 <- gamm(log_density ~ s(day, k = 10) + temp*nut, random = list(all_rep = ~1), data = eup)
summary(e.dens8$gam)
# Model Comparison
AIC(e.dens1$lme, e.dens2$lme,e.dens3$lme, e.dens4$lme,e.dens5$lme, e.dens6$lme,e.dens7$lme,e.dens8$lme)2.4 Tetrahymena Area
Same model as above, but for Tetrahymena phenotypic data.Model 7 with temperature additive and nutrient interactive effects (nutrients within smooth term) has the lowest AIC.
# 1. Temp additive effect
t.area1 <- gamm(area ~ s(day, k = 7) + temp, random = list(all_rep = ~1), data = tetraSumm,correlation = corARMA(form = ~ day|temp/all_rep, p = 3))
summary(t.area1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 7) + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 766.854 9.275 82.676 < 2e-16 ***
## temp25C -37.762 13.643 -2.768 0.00626 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.564 5.564 175.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.725
## Scale est. = 6596.5 n = 180## [1] "factor"# 2. Nut additive effect
t.area2 <- gamm(area ~ s(day, k = 7) + nut,random = list(all_rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|nut/all_rep, p = 3))
summary(t.area2$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 7) + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 773.986 8.401 92.134 < 2e-16 ***
## nutHalf -48.153 11.800 -4.081 6.85e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.616 5.616 167.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.74
## Scale est. = 6545.7 n = 180# 3. Temp and Nut additive effects
t.area3 <- gamm(area ~ s(day, k = 7) + nut + temp,random = list(all_rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|all_rep, p = 3))
summary(t.area3$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 7) + nut + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 792.262 7.824 101.255 < 2e-16 ***
## nutHalf -49.001 9.260 -5.292 3.67e-07 ***
## temp25C -38.795 9.796 -3.960 0.00011 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.601 5.601 164.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.754
## Scale est. = 6512.6 n = 180# 4. Temp interactive effect ONLY
t.area4 <- gamm(area ~ s(day, by = temp, k = 7) + temp, random = list(all_rep = ~1), data = tetraSumm, correlation = corARMA(form = ~ day|temp/all_rep, p = 3))
summary(t.area4$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = temp, k = 7) + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 769.705 9.541 80.672 <2e-16 ***
## temp25C -48.219 22.899 -2.106 0.0367 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):temp22C 5.504 5.504 178.92 <2e-16 ***
## s(day):temp25C 3.778 3.778 98.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.726
## Scale est. = 6595 n = 180# 5. Nut interactive effect ONLY
t.area5 <- gamm(area ~ s(day, by = nut, k = 7) + nut, random = list(all_rep = ~1), data = tetraSumm, correlation = corARMA(form = ~ day|nut/all_rep, p = 3))
summary(t.area5$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = nut, k = 7) + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 777.608 9.138 85.094 < 2e-16 ***
## nutHalf -55.164 13.089 -4.214 4.08e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):nutFull 4.954 4.954 102.58 <2e-16 ***
## s(day):nutHalf 5.205 5.205 68.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.766
## Scale est. = 5534.4 n = 180# 6. Temp interactive + nut additive
t.area6 <- gamm(area ~ s(day, by = temp, k = 7) + nut + temp, random = list(all_rep = ~1), data = tetraSumm, correlation = corARMA(form = ~ day|all_rep, p = 3))
summary(t.area6$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = temp, k = 7) + nut + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 794.046 8.046 98.688 < 2e-16 ***
## nutHalf -48.182 9.290 -5.186 6.12e-07 ***
## temp25C -45.377 20.821 -2.179 0.0307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):temp22C 5.536 5.536 168.60 <2e-16 ***
## s(day):temp25C 3.808 3.808 90.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.757
## Scale est. = 6493 n = 180# 7. Nut interactive + temp additive
t.area7 <- gamm(area ~ s(day, by = nut, k = 7) + temp + nut, random = list(all_rep = ~1), data = tetraSumm, correlation = corARMA(form = ~ day|all_rep, p =3))
summary(t.area7$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = nut, k = 7) + temp + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 798.259 8.251 96.751 < 2e-16 ***
## temp25C -42.354 10.314 -4.106 6.28e-05 ***
## nutHalf -56.549 9.929 -5.696 5.44e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):nutFull 4.980 4.980 111.76 <2e-16 ***
## s(day):nutHalf 5.304 5.304 77.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.783
## Scale est. = 5394 n = 180# 8. Nut * temp
t.area8 <- gamm(area ~ s(day, k = 7) + temp * nut, random = list(all_rep = ~1), data = tetraSumm, correlation = corARMA(form = ~ day|all_rep, p = 3))
summary(t.area8$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 7) + temp * nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 786.566 8.305 94.712 < 2e-16 ***
## temp25C -23.485 12.864 -1.826 0.06965 .
## nutHalf -36.548 11.660 -3.135 0.00203 **
## temp25C:nutHalf -32.309 17.733 -1.822 0.07021 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.599 5.599 153.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.757
## Scale est. = 6486 n = 180# model comparison
AIC(t.area1$lme, t.area2$lme,t.area3$lme, t.area4$lme,t.area5$lme, t.area6$lme,t.area7$lme, t.area8$lme)## df AIC
## t.area1$lme 9 2141.414
## t.area2$lme 9 2135.699
## t.area3$lme 10 2125.098
## t.area4$lme 11 2152.360
## t.area5$lme 11 2131.108
## t.area6$lme 12 2136.618
## t.area7$lme 12 2119.911
## t.area8$lme 11 2124.0722.5 Euplotes Area
Same model as above, but for Euplotes phenotypic data. Model 6 with nutrient additive and temperature interactive effects (temperature within smooth term) is the most parsimonious model.
# 1. Temp additive effect
eup.area1 <- gamm(area ~ s(day, k = 9) + temp, random = list(all_rep = ~1), data = eupSummtrait,correlation = corARMA(form = ~ day|temp/all_rep, p = 1))
summary(eup.area1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 9) + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7773.5 102.6 75.762 < 2e-16 ***
## temp25C -566.0 143.9 -3.935 0.000107 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.353 5.353 40.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.469
## Scale est. = 1.318e+06 n = 267# 2. Nut additive effect
eup.area2 <- gamm(area ~ s(day, k = 9) + nut,random = list(all_rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|nut/all_rep, p = 1))
summary(eup.area2$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 9) + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7612.9 116.2 65.492 <2e-16 ***
## nutHalf -251.5 163.7 -1.536 0.126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.323 5.323 38.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.442
## Scale est. = 1.3584e+06 n = 267# 3. Temp and Nut additive effects
eup.area3 <- gamm(area ~ s(day, k = 9) + nut + temp,random = list(all_rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(eup.area3$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 9) + nut + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7903.3 125.1 63.181 < 2e-16 ***
## nutHalf -252.9 142.6 -1.774 0.0772 .
## temp25C -569.7 142.6 -3.996 8.42e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.372 5.372 40.83 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.474
## Scale est. = 1.3016e+06 n = 267# 4. Temp interactive effect ONLY
eup.area4 <- gamm(area ~ s(day, by = temp, k = 9) + temp, random = list(all_rep = ~1), data = eupSummtrait, correlation = corARMA(form = ~ day|temp/all_rep, p = 1))
summary(eup.area4$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = temp, k = 9) + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7768.00 94.05 82.593 < 2e-16 ***
## temp25C -569.84 131.70 -4.327 2.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):temp22C 4.701 4.701 31.19 <2e-16 ***
## s(day):temp25C 5.291 5.291 24.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.515
## Scale est. = 1.1995e+06 n = 267# 5. Nut interactive effect ONLY
eup.area5 <- gamm(area ~ s(day, by = nut, k = 9) + nut, random = list(all_rep = ~1), data = eupSummtrait, correlation = corARMA(form = ~ day|nut/all_rep, p = 1))
summary(eup.area5$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = nut, k = 9) + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7606.7 116.3 65.411 <2e-16 ***
## nutHalf -244.8 163.8 -1.495 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):nutFull 4.418 4.418 24.82 <2e-16 ***
## s(day):nutHalf 5.216 5.216 24.86 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.463
## Scale est. = 1.2829e+06 n = 267# 6. Temp interactive + nut additive
eup.area6 <- gamm(area ~ s(day, by = temp, k = 9) + nut + temp, random = list(all_rep = ~1), data = eupSummtrait, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(eup.area6$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = temp, k = 9) + nut + temp
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7890.2 115.0 68.616 < 2e-16 ***
## nutHalf -237.7 130.9 -1.816 0.0705 .
## temp25C -573.2 130.7 -4.387 1.69e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):temp22C 4.702 4.702 31.26 <2e-16 ***
## s(day):temp25C 5.331 5.331 24.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.52
## Scale est. = 1.1839e+06 n = 267# 7. Nut interactive + temp additive
eup.area7 <- gamm(area ~ s(day, by = nut, k = 9) + temp + nut, random = list(all_rep = ~1), data = eupSummtrait, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(eup.area7$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, by = nut, k = 9) + temp + nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7901.6 116.2 68.001 < 2e-16 ***
## temp25C -583.0 132.2 -4.409 1.53e-05 ***
## nutHalf -243.8 132.0 -1.847 0.0659 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day):nutFull 4.438 4.438 25.83 <2e-16 ***
## s(day):nutHalf 5.332 5.332 27.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.496
## Scale est. = 1.2416e+06 n = 267# 8. nut * temp outside of smooth term
eup.area8 <- gamm(area ~ s(day, k = 9) + temp *nut, random = list(all_rep = ~1), data = eupSummtrait, correlation = corARMA(form = ~ day|all_rep, p = 1))
summary(eup.area8$gam) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## area ~ s(day, k = 9) + temp * nut
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7916.28 145.74 54.317 < 2e-16 ***
## temp25C -594.92 203.32 -2.926 0.00374 **
## nutHalf -278.24 203.91 -1.364 0.17361
## temp25C:nutHalf 49.82 286.00 0.174 0.86184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(day) 5.369 5.369 40.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.472
## Scale est. = 1.3015e+06 n = 267# model comparison
AIC(eup.area1$lme, eup.area2$lme,eup.area3$lme, eup.area4$lme,eup.area5$lme, eup.area6$lme,eup.area7$lme, eup.area8$lme)## df AIC
## eup.area1$lme 7 4551.526
## eup.area2$lme 7 4562.876
## eup.area3$lme 8 4550.373
## eup.area4$lme 9 4544.981
## eup.area5$lme 9 4569.363
## eup.area6$lme 10 4543.674
## eup.area7$lme 10 4554.992
## eup.area8$lme 9 4552.3423. Appendix I Table S3. Linear Models for Ecological Dynamics Descriptor
3.1 Bacteria OD600
The following are ecological descriptor calculated for each species.
## New data frame
bac_desc <- data.frame("tre" = rep(c("Full22C", "Full25C", "Half22C", "Half25C"), each = 6))
bac_desc<- bac_desc %>%
mutate(temp = as.factor(rep(c(22,25), each=6, times=2)),
nut = factor(rep(c("Full","Half"), each=12), levels = c("Half", "Full")),
rep = rep(1:6, times = 4))
## initial growth rate of Bacteria from day0 to day 1
bac_day0 <- bac[which(bac$day == 0), ]
bac_day1 <- bac[which(bac$day == 1), ]
bac_desc$init_growth <- (bac_day1$OD-bac_day0$OD)/8
## max abundance: three days with the top three population density
bac_max1<- teb.mean%>% group_by(tre,treatNum) %>% select(tre:day,Bacteria) %>% na.omit() %>% arrange(desc(Bacteria), .by_group = TRUE)%>% slice_max(Bacteria, n=3)
bm1 <- bac %>% filter(tre == "Full22C", day %in% bac_max1$day[which(bac_max1$tre=="Full22C")])
bm2 <- bac %>% filter(tre == "Full25C", day %in% bac_max1$day[which(bac_max1$tre=="Full25C")])
bm3 <- bac %>% filter(tre == "Half22C", day %in% bac_max1$day[which(bac_max1$tre=="Half22C")])
bm4 <- bac %>% filter(tre == "Half25C", day %in% bac_max1$day[which(bac_max1$tre=="Half25C")])
bac_max <- rbind(bm1, bm2, bm3,bm4) %>% select(tre:temp, OD, -Date)
bac_max$nut <- factor(bac_max$nut, levels = c("Half", "Full"))
# CV
bac_desc<- bac%>% group_by(tre,rep) %>%na.omit()%>% summarize(sd = sd(OD), meanOD= mean(OD), CV = sd/mean) %>% full_join(bac_desc)
# mean
bac.mean_desc <- bac_desc %>%
group_by(tre, nut, temp) %>%
na.omit %>%
summarise(ig = mean(init_growth),
CV.m=mean(CV))
bac_max_mean <- bac_max %>% group_by(nut, temp,tre) %>% na.omit() %>% summarise(ODmax = mean(OD),ODse = sd(OD))3.1.1 Bacteria Initial growth Linear Models
b.ig <- lm(init_growth ~ temp , data=bac_desc);summary(b.ig)
b.ig1 <- lm(init_growth ~ nut, data=bac_desc);summary(b.ig1)
b.ig2 <- lm(init_growth ~ temp + nut, data=bac_desc);summary(b.ig2) # no significant effect from neither temp nor nut
b.ig3 <- lm(init_growth ~ temp * nut, data=bac_desc);summary(b.ig3)
#anova(b.ig, b.ig1,b.ig2, b.ig3)
AIC(b.ig, b.ig1,b.ig2, b.ig3)3.1.2 Bacteria Maximum abundance Linear Models
b.max <- lm(OD ~ temp, bac_max);summary(b.max) # significant temperature effects
b.max1 <- lm(OD ~ nut, bac_max);summary(b.max1) # Significant nutrient effect
b.max2 <- lm(OD ~ nut+temp, bac_max);summary(b.max2) # significant nutrient and temp effect
b.max3 <- lm(OD ~ nut*temp, bac_max);summary(b.max3) # no interaction
anova(b.max, b.max2,b.max1,b.max3)
AIC(b.max, b.max1,b.max2,b.max3)3.1.3 Bacteria Coefficient of variation (CV) Linear Models
# coefficient of variation
b.cv <- lm(CV ~ temp, data=bac_desc);summary(b.cv) # temp not significant
b.cv1 <- lm(CV ~ nut, data=bac_desc);summary(b.cv1) # nut significant
b.cv2<- lm(CV ~ nut+temp, data=bac_desc);summary(b.cv2) # nut significant
b.cv3 <- lm(CV ~ nut*temp, data=bac_desc);summary(b.cv3) # nut significant
anova(b.cv, b.cv1,b.cv2, b.cv3) # two model significantly different. Interactive model has higher R2
AIC(b.cv, b.cv1,b.cv2, b.cv3)3.2 Tetrahymena
# new data frame
tetra_desc <- data.frame("tre" = rep(c("Full22C", "Full25C", "Half22C", "Half25C"), each = 6))
tetra_desc<-tetra_desc %>%
mutate(temp = as.factor(rep(c(22,25), each=6, times=2)),
nut = factor(rep(c("Full","Half"), each=12), levels = c("Half", "Full")),
rep = rep(1:6, times = 4))
## initial growth rate of tetrahymena from day0 to day 1
tetra_day0 <- tetra%>% filter(day==0) %>% select(density)
tetra_day1 <- tetra%>% filter(day==1) %>% select(density)
tetra_desc$tetra.ig <- log(tetra_day1$density)-log(tetra_day0$density)
## max abundance
tetra_max1 <- teb.mean %>% group_by(tre,treatNum) %>% select(tre:day,Tetrahymena) %>% summarize(max = max(Tetrahymena))
tetra_max <- teb.mean %>% filter(Tetrahymena %in% tetra_max1$max)
tetra_desc <- tetra%>% filter(tre %in% tetra_max$tre & day %in% tetra_max$day)%>% select(tre,rep,density) %>% rename(max = density)%>% left_join(tetra_desc)
# CV
tetra_desc <- tetra%>% group_by(tre,rep) %>% summarize(sd = sd(density), mean= mean(density), CV = sd/mean)%>% full_join(tetra_desc)
## time collapse
tetra_desc <- tetra %>% group_by(tre,rep) %>% filter(density == 0) %>% summarize(tc = min(day)) %>% full_join(tetra_desc)
# mean collapse time per treatment
tetra_desc %>% na.omit%>% group_by(tre) %>% summarize(mean_tc = mean(tc))## # A tibble: 4 × 2
## tre mean_tc
## <chr> <dbl>
## 1 Full22C 11.5
## 2 Full25C 8.5
## 3 Half22C 10.6
## 4 Half25C 8.17# mean of the descriptors
t.mean_desc <- tetra_desc %>%
group_by(tre, nut, temp) %>%
summarise(ig = mean(tetra.ig),
max.m = mean(max),
CV.m = mean(CV),
tc.m = mean(tc)) 3.2.1 Tetrahymena Initial growth Linear Models
# initial growth
t.ig <- lm(tetra.ig ~ temp, data=tetra_desc); summary(t.ig) #temp significant
t.ig1 <- lm(tetra.ig ~ nut, data=tetra_desc); summary(t.ig1) # nut NOT significant
t.ig2 <- lm(tetra.ig ~ temp + nut, data=tetra_desc); summary(t.ig2) # significant temp
t.ig3 <- lm(tetra.ig ~ temp * nut, data=tetra_desc);summary(t.ig3) # significant temp no interaction
AIC(t.ig,t.ig1, t.ig2,t.ig3)3.2.2 Tetrahymena Maximum abundance Linear Models
# max abundance
t.max <- lm(max ~ temp, data=tetra_desc); summary(t.max)
t.max1 <- lm(max ~ nut , data=tetra_desc); summary(t.max1) # significant nutrient
t.max2 <- lm(max ~ temp + nut, data=tetra_desc); summary(t.max2) # significant nutrient, NO temp
t.max3 <- lm(max ~ temp * nut, data=tetra_desc); summary(t.max3) # significant temp and nutrient and interactive effect. Best model
anova(t.max2, t.max1) # model significantly different
AIC(t.max,t.max1, t.max2,t.max3)3.2.3 Tetrahymena Coefficient of variation (CV) Linear Models
#Coefficient of varitaion
t.cv <- lm(CV ~ temp, data=tetra_desc);summary(t.cv)
t.cv1 <- lm(CV ~ nut, data=tetra_desc);summary(t.cv1)
t.cv2 <- lm(CV ~ nut+temp, data=tetra_desc);summary(t.cv2) # Temp significant. nut p = 0.08
t.cv3 <- lm(CV ~ nut*temp, data=tetra_desc);summary(t.cv3) # temp, nut and interaction significant, 0.08
anova(t.cv1, t.cv2) # two model significantly different. interaction model better
AIC(t.cv, t.cv1,t.cv2, t.cv3)3.2.4 Tetrahymena Day of population collapse Linear Models
# time extinction/ day of collapse
tExt <- lm(tc ~ temp, data=tetra_desc); summary(tExt) # significant Temperature effects
tExt1 <- lm(tc~ nut, data=tetra_desc);summary(tExt1)
tExt2 <- lm(tc~ temp + nut, data=tetra_desc);summary(tExt2)
tExt3 <- lm(tc~ temp * nut, data=tetra_desc);summary(tExt3)
anova(tExt1, tExt2)
AIC(tExt, tExt1, tExt2, tExt3)3.3 Euplotes
3.3.1 Euplotes Initial growth Linear Models
# Initial Growth
e.ig <- lm(init_growth ~ temp, data=eup_desc);summary(e.ig)
e.ig1 <- lm(init_growth ~ nut, data=eup_desc);summary(e.ig1) # no nut
e.ig2 <- lm(init_growth ~ temp + nut, data=eup_desc);summary(e.ig2) # significant temp
e.ig3 <- lm(init_growth ~ temp * nut, data=eup_desc);summary(e.ig3) # significant temp, sig nut, and interactive!
anova(e.ig, e.ig1,e.ig2, e.ig3)
AIC(e.ig, e.ig1,e.ig2, e.ig3)3.3.2 Euplotes Maximum abundance Linear Models
# Maximum Abundance
e.max <- lm(max ~ temp, data=eup_desc);summary(e.max)
e.max1 <- lm(max ~ nut, data=eup_desc);summary(e.max1)
e.max2 <- lm(max ~ nut+temp, data=eup_desc);summary(e.max2) # significant temp and nut additive
e.max3 <- lm(max ~ nut*temp, data=eup_desc);summary(e.max3) # no interactive
anova(e.max2, e.max3)
AIC(e.max, e.max1, e.max2, e.max3)3.3.3 Euplotes Coefficient of variation (CV) Linear Models
# coefficient of variation
e.cv <- lm(CV ~ temp, data=eup_desc);summary(e.cv)
e.cv1 <- lm(CV ~ nut, data=eup_desc);summary(e.cv1)
e.cv2 <- lm(CV ~ nut+temp, data=eup_desc);summary(e.cv2) # nut and temp both significant
e.cv3 <- lm(CV ~ nut*temp, data=eup_desc);summary(e.cv3) # nut, temp and interactive all significant R2 = 0.7545
anova(e.cv1, e.cv2) # two model significantly different. interaction model better
AIC(e.cv,e.cv1, e.cv2,e.cv3)3.3.4 Euplotes Day of population peak Linear Models
# day peaked
e.dp <- lm(daypeak ~ temp, data=eup_desc);summary(e.dp)
e.dp1 <- lm(daypeak ~ nut, data=eup_desc);summary(e.dp1)
e.dp2 <- lm(daypeak ~ nut+temp, data=eup_desc);summary(e.dp2) # significant temp, no nut
e.dp3 <- lm(daypeak ~ nut*temp, data=eup_desc);summary(e.dp3) # significant temp, no interavtive nor nut
anova(e.dp1, e.dp2, e.dp3)
AIC(e.dp,e.dp1, e.dp2, e.dp3)4. Appendix II Figure S1-14 CCM convergence plot for each species pair
We started by generating GAMM predictions from the selected model as the time series used in CCM test.
# function to predict tetrahymena, euplotes, and bacteria density from previous selected models
em0 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = eup)
tm1 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 1))
bm2 <- gamm(OD ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 2))
pred_gam_TEB <- function(treatment,day, length){
df <- data.frame(day=seq(0, day,length.out=length), tre=rep(treatment,length))
eup_pred <- exp(predict.gam(em0$gam,df,type = "response", se.fit =T)$fit)-1
tetra_pred <- exp(predict.gam(tm1$gam,df,type = "response", se.fit =T)$fit)-1
bac_pred <- predict.gam(bm2$gam,df,type = "response", se.fit =T)$fit
pred <- as.data.frame(cbind(df[, 1], eup_pred, tetra_pred, bac_pred))
colnames(pred) <- c("day","eup_pred", "tetra_pred", "bac_pred" )
pred
}
F22gam <- pred_gam_TEB("Full22C",15, 800)
F25gam <- pred_gam_TEB("Full25C",15, 800)
H22gam <- pred_gam_TEB("Half22C",15, 800)
H25gam <- pred_gam_TEB("Half25C",15, 800)
# cropped time series for ecological dynamics
F22gam_crop <- pred_gam_TEB("Full22C",9, 500)
F25gam_crop <- pred_gam_TEB("Full25C",9, 500)
H22gam_crop <- pred_gam_TEB("Half22C",9, 500)
H25gam_crop <- pred_gam_TEB("Half25C",9, 500)# selected gamm model fpr tetrahymena and Euplotes trait
tt3 <- gamm(area ~ s(day, by = tre, k =7) + tre,random = list(rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|tre/rep, p = 3))
et1 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSummtrait, correlation = corARMA(form = ~ day|tre/rep, p = 1))
# predict Gamm function for trait
pred_gam_TEB_trait <- function(day, tModel, eModel, treatment, length, densityData){
df <- data.frame(day=seq(0,day,length.out=length), tre=rep(treatment,length))
eup_trait <- predict.gam(eModel$gam,df,type = "response", se.fit =T)$fit
tetra_trait <- predict.gam(tModel$gam,df,type = "response", se.fit =T)$fit
pred <- as.data.frame(cbind(densityData,eup_trait,tetra_trait))
}
# gam fit that goes into CCM
F22gam.area <- pred_gam_TEB_trait(15,tt3,et1,"Full22C", 800, F22gam)
F25gam.area <- pred_gam_TEB_trait(15,tt3,et1,"Full25C", 800, F25gam)
H22gam.area <- pred_gam_TEB_trait(15,tt3,et1,"Half22C", 800, H22gam)
H25gam.area <- pred_gam_TEB_trait(15,tt3,et1,"Half25C", 800, H25gam)
# Cropped time series for phenotypic dynamics, with cropped predicted ecological dynamics merged
F22gam.area_crop <- pred_gam_TEB_trait(9,tt3,et1,"Full22C", 500, F22gam_crop)
F25gam.area_crop <- pred_gam_TEB_trait(9,tt3,et1,"Full25C", 500, F25gam_crop)
H22gam.area_crop <- pred_gam_TEB_trait(9,tt3,et1,"Half22C", 500, H22gam_crop)
H25gam.area_crop <- pred_gam_TEB_trait(9,tt3,et1,"Half25C", 500, H25gam_crop)There are two seperate functions for calculating the CCM predictability. All CCM pairs with tetrahymena trait datas, the model fits only include gamm model prediction from day 0 to day 9, since all Tetrahymena populations collaplse around day 9 based on gamm prediction.
# this function only calculate the CCM analysis involve tetrahymena phenotypic dynamics. Here we cropped all the time series involving in these CCM.
CCM_TEB_trait_crop <- function(data, Em, libSize, nSamples){
E_Xmap_Tt <- CCM(dataFrame=data, E=Em, columns = "eup_pred" , target = "tetra_trait", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
Et_Xmap_Tt <- CCM(dataFrame=data, E=Em, columns = "eup_trait" , target = "tetra_trait", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
T_Xmap_Tt <- CCM(dataFrame=data, E=Em, columns = "tetra_pred" , target = "tetra_trait", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
Tt_Xmap_B <- CCM(dataFrame=data, E=Em, columns = "tetra_trait" , target = "bac_pred", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
a <- list( "E_Xmap_Tt" = E_Xmap_Tt[,-3], "Tt_Xmap_E" = E_Xmap_Tt[,-2], "Et_Xmap_Tt" = Et_Xmap_Tt[,-3], "Tt_Xmap_Et"= Et_Xmap_Tt[,-2], "T_Xmap_Tt" = T_Xmap_Tt[,-3], "Tt_Xmap_T" = T_Xmap_Tt[,-2], "Tt_Xmap_B" = Tt_Xmap_B[,-3], "B_Xmap_Tt" = Tt_Xmap_B[,-2])
}
# this function calculate all the rest CCM pairs that do not involve tetra trait data
# dateset from model prediction used in this function are different length from the ones used in the previous function
CCM_TEB_trait_full <- function(data, Em, libSize, nSamples){
E_Xmap_T <- CCM(dataFrame=data, E=Em, columns ="eup_pred" , target = "tetra_pred", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
T_Xmap_B <- CCM(dataFrame=data, E=Em, columns ="tetra_pred" , target = "bac_pred", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
E_Xmap_B <- CCM(dataFrame=data, E=Em, columns ="eup_pred" , target = "bac_pred",libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
Et_Xmap_T <- CCM(dataFrame=data, E=Em, columns ="eup_trait" ,target = "tetra_pred",libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
E_Xmap_Et <- CCM(dataFrame=data, E=Em, columns = "eup_pred", target = "eup_trait", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
Et_Xmap_B <- CCM(dataFrame=data, E=Em, columns = "eup_trait" , target = "bac_pred", libSizes = libSize, sample = nSamples, random = TRUE, replacement = TRUE)
a <- list("E_Xmap_T" = E_Xmap_T[,-3],"T_Xmap_E" = E_Xmap_T[,-2], "T_Xmap_B" = T_Xmap_B[,-3], "B_Xmap_T" = T_Xmap_B[, -2], "E_Xmap_B" = E_Xmap_B[,-3], "B_Xmap_E"= E_Xmap_B[,-2],"E_Xmap_Et" = E_Xmap_Et[,-3], "Et_Xmap_E" = E_Xmap_Et[,-2], "Et_Xmap_T" = Et_Xmap_T[,-3], "T_Xmap_Et" = Et_Xmap_T[,-2], "Et_Xmap_B" = Et_Xmap_B[,-3],"B_Xmap_Et" = Et_Xmap_B[,-2])
}Here are mean CCM predictability calcualted from 100 replication.
# calculating CCM involving tetrahymena trait dynamics, using cropped time series
F22CCM.area_crop <-CCM_TEB_trait_crop(F22gam.area_crop, 1, seq(10, 80, by = 10),100)
F25CCM.area_crop <-CCM_TEB_trait_crop(F25gam.area_crop, 1, seq(10, 80, by = 10),100)
H22CCM.area_crop <-CCM_TEB_trait_crop(H22gam.area_crop, 1, seq(10, 80, by = 10),100)
H25CCM.area_crop <-CCM_TEB_trait_crop(H25gam.area_crop, 1, seq(10, 80, by = 10),100)
# calculate the rest CCM using full time series
F22CCM.area <-CCM_TEB_trait_full(F22gam.area, 1, seq(10, 80, by = 10),100)
F25CCM.area <-CCM_TEB_trait_full(F25gam.area, 1, seq(10, 80, by = 10),100)
H22CCM.area <-CCM_TEB_trait_full(H22gam.area, 1, seq(10, 80, by = 10),100)
H25CCM.area <-CCM_TEB_trait_full(H25gam.area, 1, seq(10, 80, by = 10),100)Here is the function for calculating the CCM predictability n times (“nSamples”) with increasing library size. This function is specific used for generating the CCM convergence plot. Due to the low counts for phenotypic data gathered near T. pyriformis population collapse time, GAMM predictions become unreliable. To better quantifying the relationship between T. pyriformis phenotypic dynamics with other species, we only used GAMM prediction from day 0-9 for T. pyriformis phenotypic dynamics. Correspondingly, we only use GAMM prediction from 0-9 for any other time series when paired up with T. pyriformis phenotypic dynamics. Therefore, there is one function desinated for all the pair of CCM test involving T. pyriformis phenotypic dynamics to used the cropped time series from day 0-9.
# Note: running E_Xmap_Tt <- ccm(..stats_only = FALSE) by itself out of the function is fine. BUT DO NOT VIEW E_Xmap_Tt. It's likely to overload R and abort the session.Check E_Xmap_Tt$CCM1_PredictStat and $CCM2_PredictStat for predictability results
# pairs with Tetrahymena phenotypic daynamics
ccmTEB_trait_info_crop <- function(data, Em, libSize, nSamples){
E_Xmap_Tt <- ccm(data, E=Em, lib_column = "eup_pred" , target_column = "tetra_trait", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
Et_Xmap_Tt <- ccm(data, E=Em, lib_column = "eup_trait" , target_column = "tetra_trait", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
T_Xmap_Tt <- ccm(data, E=Em, lib_column = "tetra_pred" , target_column = "tetra_trait", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
Tt_Xmap_B <- ccm(data, E=Em, lib_column = "tetra_trait" , target_column = "bac_pred", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
a <- list("E_Xmap_Tt" = E_Xmap_Tt$CCM1_PredictStat, "Tt_Xmap_E" = E_Xmap_Tt$CCM2_PredictStat, "Et_Xmap_Tt" = Et_Xmap_Tt$CCM1_PredictStat, "Tt_Xmap_Et"= Et_Xmap_Tt$CCM2_PredictStat, "T_Xmap_Tt" = T_Xmap_Tt$CCM1_PredictStat, "Tt_Xmap_T" = T_Xmap_Tt$CCM2_PredictStat, "Tt_Xmap_B" = Tt_Xmap_B$CCM1_PredictStat, "B_Xmap_Tt" = Tt_Xmap_B$CCM2_PredictStat)
}
# pairs that do not involve Tetrahymena phenotypic daynamics
ccmTEB_trait_info <- function(data, Em, libSize, nSamples){
E_Xmap_T <- ccm(data, E=Em, lib_column ="eup_pred" , target_column = "tetra_pred", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
T_Xmap_B <- ccm(data, E=Em, lib_column ="tetra_pred" , target_column = "bac_pred", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
E_Xmap_B <- ccm(data, E=Em, lib_column ="eup_pred" , target_column = "bac_pred",lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
Et_Xmap_T <- ccm(data, E=Em, lib_column ="eup_trait" ,target_column = "tetra_pred",lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
E_Xmap_Et <- ccm(data, E=Em, lib_column = "eup_pred", target_column = "eup_trait", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
Et_Xmap_B <- ccm(data, E=Em, lib_column = "eup_trait" , target_column = "bac_pred", lib_sizes = libSize, num_samples = nSamples, random_libs = TRUE, replace = TRUE, stats_only = FALSE)
a <- list("E_Xmap_T" = E_Xmap_T$CCM1_PredictStat, "T_Xmap_E" = E_Xmap_T$CCM2_PredictStat, "T_Xmap_B" = T_Xmap_B$CCM1_PredictStat, "B_Xmap_T" = T_Xmap_B$CCM2_PredictStat, "E_Xmap_B" = E_Xmap_B$CCM1_PredictStat, "B_Xmap_E"= E_Xmap_B$CCM2_PredictStat,"E_Xmap_Et" = E_Xmap_Et$CCM1_PredictStat, "Et_Xmap_E" = E_Xmap_Et$CCM2_PredictStat, "Et_Xmap_T" = Et_Xmap_T$CCM1_PredictStat, "T_Xmap_Et" = Et_Xmap_T$CCM2_PredictStat, "Et_Xmap_B" = Et_Xmap_B$CCM1_PredictStat,"B_Xmap_Et" = Et_Xmap_B$CCM2_PredictStat)
}
F22ccmtrait_info1 <- ccmTEB_trait_info_crop(F22gam.area_crop, 1, seq(10, 80, by = 10),100) # changing the library size to 200 doesn't qualitatively change the predictability
F25ccmtrait_info1 <- ccmTEB_trait_info_crop(F25gam.area_crop, 1, seq(10, 80, by = 10),100)
H22ccmtrait_info1 <- ccmTEB_trait_info_crop(H22gam.area_crop, 1, seq(10, 80, by = 10),100)
H25ccmtrait_info1 <- ccmTEB_trait_info_crop(H25gam.area_crop, 1, seq(10, 80, by = 10),100)
F22ccmtrait_info2 <- ccmTEB_trait_info(F22gam.area, 1, seq(10, 80, by = 10),100) # changing the library size to 200 doesn't qualitatively change the predictability
F25ccmtrait_info2 <- ccmTEB_trait_info(F25gam.area, 1, seq(10, 80, by = 10),100)
H22ccmtrait_info2 <- ccmTEB_trait_info(H22gam.area, 1, seq(10, 80, by = 10),100)
H25ccmtrait_info2 <- ccmTEB_trait_info(H25gam.area, 1, seq(10, 80, by = 10),100)4.1 Eco-Eco Covergence Plot
# Half 22C
par(mfrow = c(3,2))
plot(rho~LibSize, H22ccmtrait_info2$E_Xmap_T, ylim = c(0,1),main = "L22 Tetra Eco->Eup Eco ", col = "#cfd3db");lines(H22CCM.area$E_Xmap_T[,1], H22CCM.area$E_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$T_Xmap_E, ylim = c(0,1), main = "L22 Eup Eco->Tetra Eco", col = "#cfd3db");lines(H22CCM.area$T_Xmap_E[,1], H22CCM.area$T_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info2$E_Xmap_B, ylim = c(0,1), main = "L22 Bac Eco-> Eup Eco", col = "#cfd3db");lines(H22CCM.area$E_Xmap_B[,1], H22CCM.area$E_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$B_Xmap_E, ylim = c(0,1), main = "L22 Eup Eco->Bac Eco", col = "#cfd3db");lines(H22CCM.area$B_Xmap_E[,1], H22CCM.area$B_Xmap_E[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info2$T_Xmap_B, ylim = c(0,1), main = "L22 Bac Eco->Tetra Eco", col = "#cfd3db");lines(H22CCM.area$T_Xmap_B[,1], H22CCM.area$T_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$B_Xmap_T, ylim = c(0,1), main = "L22 Tetra Eco->Bac Eco", col = "#cfd3db");lines(H22CCM.area$B_Xmap_T[,1], H22CCM.area$B_Xmap_T[,2], col = "#A997DF", lwd = 2, lty = 2)
# Full 22C
par(mfrow = c(3,2))
plot(rho~LibSize, F22ccmtrait_info2$E_Xmap_T, ylim = c(0,1),main = "F22 Tetra Eco->Eup Eco ", col = "#cfd3db");lines(F22CCM.area$E_Xmap_T[,1], F22CCM.area$E_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$T_Xmap_E, ylim = c(0,1), main = "F22 Eup Eco->Tetra Eco", col = "#cfd3db");lines(F22CCM.area$T_Xmap_E[,1], F22CCM.area$T_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info2$E_Xmap_B, ylim = c(0,1), main = "F22 Bac Eco-> Eup Eco", col = "#cfd3db");lines(F22CCM.area$E_Xmap_B[,1], F22CCM.area$E_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$B_Xmap_E, ylim = c(0,1), main = "F22 Eup Eco->Bac Eco", col = "#cfd3db");lines(F22CCM.area$B_Xmap_E[,1], F22CCM.area$B_Xmap_E[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info2$T_Xmap_B, ylim = c(0,1), main = "F22 Bac Eco->Tetra Eco", col = "#cfd3db");lines(F22CCM.area$T_Xmap_B[,1], F22CCM.area$T_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$B_Xmap_T, ylim = c(0,1), main = "F22 Tetra Eco->Bac Eco", col = "#cfd3db");lines(F22CCM.area$B_Xmap_T[,1], F22CCM.area$B_Xmap_T[,2], col = "#A997DF", lwd = 2, lty = 2)
# Half 25C
par(mfrow = c(3,2))
plot(rho~LibSize, H25ccmtrait_info2$E_Xmap_T, ylim = c(0,1),main = "H25 Tetra Eco->Eup Eco ", col = "#cfd3db");lines(H25CCM.area$E_Xmap_T[,1], H25CCM.area$E_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$T_Xmap_E, ylim = c(0,1), main = "H25 Eup Eco->Tetra Eco", col = "#cfd3db");lines(H25CCM.area$T_Xmap_E[,1], H25CCM.area$T_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info2$E_Xmap_B, ylim = c(0,1), main = "H25 Bac Eco-> Eup Eco", col = "#cfd3db");lines(H25CCM.area$E_Xmap_B[,1], H25CCM.area$E_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$B_Xmap_E, ylim = c(0,1), main = "H25 Eup Eco->Bac Eco", col = "#cfd3db");lines(H25CCM.area$B_Xmap_E[,1], H25CCM.area$B_Xmap_E[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info2$T_Xmap_B, ylim = c(0,1), main = "H25 Bac Eco->Tetra Eco", col = "#cfd3db");lines(H25CCM.area$T_Xmap_B[,1], H25CCM.area$T_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$B_Xmap_T, ylim = c(0,1), main = "H25 Tetra Eco->Bac Eco", col = "#cfd3db");lines(H25CCM.area$B_Xmap_T[,1], H25CCM.area$B_Xmap_T[,2], col = "#A997DF", lwd = 2, lty = 2)
# Full 25C
par(mfrow = c(3,2))
plot(rho~LibSize, F25ccmtrait_info2$E_Xmap_T, ylim = c(0,1),main = "F25 Tetra Eco->Eup Eco ", col = "#cfd3db");lines(F25CCM.area$E_Xmap_T[,1], F25CCM.area$E_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$T_Xmap_E, ylim = c(0,1), main = "F25 Eup Eco->Tetra Eco", col = "#cfd3db");lines(F25CCM.area$T_Xmap_E[,1], F25CCM.area$T_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info2$E_Xmap_B, ylim = c(0,1), main = "F25 Bac Eco-> Eup Eco", col = "#cfd3db");lines(F25CCM.area$E_Xmap_B[,1], F25CCM.area$E_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$B_Xmap_E, ylim = c(0,1), main = "F25 Eup Eco->Bac Eco", col = "#cfd3db");lines(F25CCM.area$B_Xmap_E[,1], F25CCM.area$B_Xmap_E[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info2$T_Xmap_B, ylim = c(0,1), main = "F25 Bac Eco->Tetra Eco", col = "#cfd3db");lines(F25CCM.area$T_Xmap_B[,1], F25CCM.area$T_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$B_Xmap_T, ylim = c(0,1), main = "F25 Tetra Eco->Bac Eco", col = "#cfd3db");lines(F25CCM.area$B_Xmap_T[,1], F25CCM.area$B_Xmap_T[,2], col = "#A997DF", lwd = 2, lty = 2)
4.2 Eco-Pheno Covergence Plot
# Half 22C
par(mfrow = c(3,2))
plot(rho~LibSize, H22ccmtrait_info1$Tt_Xmap_E, ylim = c(0,1),main = "H22 Eup Eco->Tetra Pheno ", col = "#cfd3db");lines(H22CCM.area_crop$Tt_Xmap_E[,1], H22CCM.area_crop$Tt_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info2$Et_Xmap_T, ylim = c(0,1), main = "H22 Tetra Eco-> Eup Pheno", col = "#cfd3db");lines(H22CCM.area$Et_Xmap_T[,1], H22CCM.area$Et_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$Et_Xmap_B, ylim = c(0,1), main = "H22 Bac Eco-> Eup Pheno", col = "#cfd3db");lines(H22CCM.area$Et_Xmap_B[,1], H22CCM.area$Et_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info1$Tt_Xmap_B, ylim = c(0,1), main = "H22 Bac Eco-> Tetra Pheno", col = "#cfd3db");lines(H22CCM.area_crop$Tt_Xmap_B[,1], H22CCM.area_crop$Tt_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$Et_Xmap_E, ylim = c(0,1), main = "H22 Eup Eco-> Eup Pheno", col = "#cfd3db");lines(H22CCM.area$Et_Xmap_E[,1], H22CCM.area$Et_Xmap_E[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, H22ccmtrait_info1$Tt_Xmap_T, ylim = c(0,1), main = "H22 Tetra Eco-> Tetra Pheno", col = "#cfd3db");lines(H22CCM.area_crop$Tt_Xmap_T[,1], H22CCM.area_crop$Tt_Xmap_T[,2], col = "#929084", lwd = 2, lty = 3)
# Full 22C
par(mfrow = c(3,2))
plot(rho~LibSize, F22ccmtrait_info1$Tt_Xmap_E, ylim = c(0,1),main = "F22 Eup Eco->Tetra Pheno ", col = "#cfd3db");lines(F22CCM.area_crop$Tt_Xmap_E[,1], F22CCM.area_crop$Tt_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info2$Et_Xmap_T, ylim = c(0,1), main = "F22 Tetra Eco-> Eup Pheno", col = "#cfd3db");lines(F22CCM.area$Et_Xmap_T[,1], F22CCM.area$Et_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$Et_Xmap_B, ylim = c(0,1), main = "F22 Bac Eco-> Eup Pheno", col = "#cfd3db");lines(F22CCM.area$Et_Xmap_B[,1], F22CCM.area$Et_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info1$Tt_Xmap_B, ylim = c(0,1), main = "F22 Bac Eco-> Tetra Pheno", col = "#cfd3db");lines(F22CCM.area_crop$Tt_Xmap_B[,1], F22CCM.area_crop$Tt_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$Et_Xmap_E, ylim = c(0,1), main = "F22 Eup Eco-> Eup Pheno", col = "#cfd3db");lines(F22CCM.area$Et_Xmap_E[,1], F22CCM.area$Et_Xmap_E[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, F22ccmtrait_info1$Tt_Xmap_T, ylim = c(0,1), main = "F22 Tetra Eco-> Tetra Pheno", col = "#cfd3db");lines(F22CCM.area_crop$Tt_Xmap_T[,1], F22CCM.area_crop$Tt_Xmap_T[,2], col = "#929084", lwd = 2, lty = 3)
# Half 25C
par(mfrow = c(3,2))
plot(rho~LibSize, H25ccmtrait_info1$Tt_Xmap_E, ylim = c(0,1),main = "H25 Eup Eco->Tetra Pheno ", col = "#cfd3db");lines(H25CCM.area_crop$Tt_Xmap_E[,1], H25CCM.area_crop$Tt_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info2$Et_Xmap_T, ylim = c(0,1), main = "H25 Tetra Eco-> Eup Pheno", col = "#cfd3db");lines(H25CCM.area$Et_Xmap_T[,1], H25CCM.area$Et_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$Et_Xmap_B, ylim = c(0,1), main = "H25 Bac Eco-> Eup Pheno", col = "#cfd3db");lines(H25CCM.area$Et_Xmap_B[,1], H25CCM.area$Et_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info1$Tt_Xmap_B, ylim = c(0,1), main = "H25 Bac Eco-> Tetra Pheno", col = "#cfd3db");lines(H25CCM.area_crop$Tt_Xmap_B[,1], H25CCM.area_crop$Tt_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$Et_Xmap_E, ylim = c(0,1), main = "H25 Eup Eco-> Eup Pheno", col = "#cfd3db");lines(H25CCM.area$Et_Xmap_E[,1], H25CCM.area$Et_Xmap_E[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, H25ccmtrait_info1$Tt_Xmap_T, ylim = c(0,1), main = "H25 Tetra Eco-> Tetra Pheno", col = "#cfd3db");lines(H25CCM.area_crop$Tt_Xmap_T[,1], H25CCM.area_crop$Tt_Xmap_T[,2], col = "#929084", lwd = 2, lty = 3)
# Full 25C
par(mfrow = c(3,2))
plot(rho~LibSize, F25ccmtrait_info1$Tt_Xmap_E, ylim = c(0,1),main = "F25 Eup Eco->Tetra Pheno ", col = "#cfd3db");lines(F25CCM.area_crop$Tt_Xmap_E[,1], F25CCM.area_crop$Tt_Xmap_E[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info2$Et_Xmap_T, ylim = c(0,1), main = "F25 Tetra Eco-> Eup Pheno", col = "#cfd3db");lines(F25CCM.area$Et_Xmap_T[,1], F25CCM.area$Et_Xmap_T[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$Et_Xmap_B, ylim = c(0,1), main = "F25 Bac Eco-> Eup Pheno", col = "#cfd3db");lines(F25CCM.area$Et_Xmap_B[,1], F25CCM.area$Et_Xmap_B[,2], col = "#add0af", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info1$Tt_Xmap_B, ylim = c(0,1), main = "F25 Bac Eco-> Tetra Pheno", col = "#cfd3db");lines(F25CCM.area_crop$Tt_Xmap_B[,1], F25CCM.area_crop$Tt_Xmap_B[,2], col = "#A997DF", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$Et_Xmap_E, ylim = c(0,1), main = "F25 Eup Eco-> Eup Pheno", col = "#cfd3db");lines(F25CCM.area$Et_Xmap_E[,1], F25CCM.area$Et_Xmap_E[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, F25ccmtrait_info1$Tt_Xmap_T, ylim = c(0,1), main = "F25 Tetra Eco-> Tetra Pheno", col = "#cfd3db");lines(F25CCM.area_crop$Tt_Xmap_T[,1], F25CCM.area_crop$Tt_Xmap_T[,2], col = "#929084", lwd = 2, lty = 3)
4.3 Pheno-Eco Covergence Plot
# Half 22C
par(mfrow = c(3,2))
plot(rho~LibSize, H22ccmtrait_info1$E_Xmap_Tt, ylim = c(0,1),main = "H22 Tetra Pheno-> Eup Eco", col = "#cfd3db");lines(H22CCM.area_crop$E_Xmap_Tt[,1], H22CCM.area_crop$E_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info2$T_Xmap_Et, ylim = c(0,1), main = "H22 Eup Pheno-> Tetra Eco", col = "#cfd3db");lines(H22CCM.area$T_Xmap_Et[,1], H22CCM.area$T_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info2$B_Xmap_Et, ylim = c(0,1), main = "H22 Eup Pheno-> Bac Eco", col = "#cfd3db");lines(H22CCM.area$B_Xmap_Et[,1], H22CCM.area$B_Xmap_Et[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info1$B_Xmap_Tt, ylim = c(0,1), main = "H22 Tetra Pheno-> Bac Eco", col = "#cfd3db");lines(H22CCM.area_crop$B_Xmap_Tt[,1], H22CCM.area_crop$B_Xmap_Tt[,2], col = "#A997DF", lwd = 2, lty = 2)
plot(rho~LibSize, H22ccmtrait_info2$E_Xmap_Et, ylim = c(0,1), main = "H22 Eup Pheno-> Eup Eco", col = "#cfd3db");lines(H22CCM.area$E_Xmap_Et[,1], H22CCM.area$E_Xmap_Et[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, H22ccmtrait_info1$T_Xmap_Tt, ylim = c(0,1), main = "H22 Tetra Pheno-> Tetra Eco", col = "#cfd3db");lines(H22CCM.area_crop$T_Xmap_Tt[,1], H22CCM.area_crop$T_Xmap_Tt[,2], col = "#929084", lwd = 2, lty = 3)
# Full 22C
par(mfrow = c(3,2))
plot(rho~LibSize, F22ccmtrait_info1$E_Xmap_Tt, ylim = c(0,1),main = "F22 Tetra Pheno-> Eup Eco", col = "#cfd3db");lines(F22CCM.area_crop$E_Xmap_Tt[,1], F22CCM.area_crop$E_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info2$T_Xmap_Et, ylim = c(0,1), main = "F22 Eup Pheno-> Tetra Eco", col = "#cfd3db");lines(F22CCM.area$T_Xmap_Et[,1], F22CCM.area$T_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info2$B_Xmap_Et, ylim = c(0,1), main = "F22 Eup Pheno-> Bac Eco", col = "#cfd3db");lines(F22CCM.area$B_Xmap_Et[,1], F22CCM.area$B_Xmap_Et[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info1$B_Xmap_Tt, ylim = c(0,1), main = "F22 Tetra Pheno-> Bac Eco", col = "#cfd3db");lines(F22CCM.area_crop$B_Xmap_Tt[,1], F22CCM.area_crop$B_Xmap_Tt[,2], col = "#A997DF", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info2$E_Xmap_Et, ylim = c(0,1), main = "F22 Eup Pheno-> Eup Eco", col = "#cfd3db");lines(F22CCM.area$E_Xmap_Et[,1], F22CCM.area$E_Xmap_Et[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, F22ccmtrait_info1$T_Xmap_Tt, ylim = c(0,1), main = "F22 Tetra Pheno-> Tetra Eco", col = "#cfd3db");lines(F22CCM.area_crop$T_Xmap_Tt[,1], F22CCM.area_crop$T_Xmap_Tt[,2], col = "#929084", lwd = 2, lty = 3)
# Half 25C
par(mfrow = c(3,2))
plot(rho~LibSize, H25ccmtrait_info1$E_Xmap_Tt, ylim = c(0,1),main = "H25 Tetra Pheno-> Eup Eco", col = "#cfd3db");lines(H25CCM.area_crop$E_Xmap_Tt[,1], H25CCM.area_crop$E_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info2$T_Xmap_Et, ylim = c(0,1), main = "H25 Eup Pheno-> Tetra Eco", col = "#cfd3db");lines(H25CCM.area$T_Xmap_Et[,1], H25CCM.area$T_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info2$B_Xmap_Et, ylim = c(0,1), main = "H25 Eup Pheno-> Bac Eco", col = "#cfd3db");lines(H25CCM.area$B_Xmap_Et[,1], H25CCM.area$B_Xmap_Et[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info1$B_Xmap_Tt, ylim = c(0,1), main = "H25 Tetra Pheno-> Bac Eco", col = "#cfd3db");lines(H25CCM.area_crop$B_Xmap_Tt[,1], H25CCM.area_crop$B_Xmap_Tt[,2], col = "#A997DF", lwd = 2, lty = 2)
plot(rho~LibSize, H25ccmtrait_info2$E_Xmap_Et, ylim = c(0,1), main = "H25 Eup Pheno-> Eup Eco", col = "#cfd3db");lines(H25CCM.area$E_Xmap_Et[,1], H25CCM.area$E_Xmap_Et[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, H25ccmtrait_info1$T_Xmap_Tt, ylim = c(0,1), main = "H25 Tetra Pheno-> Tetra Eco", col = "#cfd3db");lines(H25CCM.area_crop$T_Xmap_Tt[,1], H25CCM.area_crop$T_Xmap_Tt[,2], col = "#929084", lwd = 2, lty = 3)
# Full 25C
par(mfrow = c(3,2))
plot(rho~LibSize, F25ccmtrait_info1$E_Xmap_Tt, ylim = c(0,1),main = "F25 Tetra Pheno-> Eup Eco", col = "#cfd3db");lines(F25CCM.area_crop$E_Xmap_Tt[,1], F25CCM.area_crop$E_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info2$T_Xmap_Et, ylim = c(0,1), main = "F25 Eup Pheno-> Tetra Eco", col = "#cfd3db");lines(F25CCM.area$T_Xmap_Et[,1], F25CCM.area$T_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info2$B_Xmap_Et, ylim = c(0,1), main = "F25 Eup Pheno-> Bac Eco", col = "#cfd3db");lines(F25CCM.area$B_Xmap_Et[,1], F25CCM.area$B_Xmap_Et[,2], col = "#add0af", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info1$B_Xmap_Tt, ylim = c(0,1), main = "F25 Tetra Pheno-> Bac Eco", col = "#cfd3db");lines(F25CCM.area_crop$B_Xmap_Tt[,1], F25CCM.area_crop$B_Xmap_Tt[,2], col = "#A997DF", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info2$E_Xmap_Et, ylim = c(0,1), main = "F25 Eup Pheno-> Eup Eco", col = "#cfd3db");lines(F25CCM.area$E_Xmap_Et[,1], F25CCM.area$E_Xmap_Et[,2], col = "#2E4052", lwd = 2, lty = 3)
plot(rho~LibSize, F25ccmtrait_info1$T_Xmap_Tt, ylim = c(0,1), main = "F25 Tetra Pheno-> Tetra Eco", col = "#cfd3db");lines(F25CCM.area_crop$T_Xmap_Tt[,1], F25CCM.area_crop$T_Xmap_Tt[,2], col = "#929084", lwd = 2, lty = 3)
4.4 Pheno-Pheno
# Half 22C and Full 22C
par(mfrow = c(2,2))
plot(rho~LibSize, H22ccmtrait_info1$Et_Xmap_Tt, ylim = c(0,1),main = "H22 Tetra Pheno-> Eup Pheno", col = "#cfd3db");lines(H22CCM.area_crop$Et_Xmap_Tt[,1], H22CCM.area_crop$Et_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H22ccmtrait_info1$Tt_Xmap_Et, ylim = c(0,1),main = "H22 Eup Pheno-> Tetra Pheno", col = "#cfd3db");lines(H22CCM.area_crop$Tt_Xmap_Et[,1], H22CCM.area_crop$Tt_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F22ccmtrait_info1$Et_Xmap_Tt, ylim = c(0,1),main = "F22 Tetra Pheno-> Eup Pheno", col = "#cfd3db");lines(F22CCM.area_crop$Et_Xmap_Tt[,1], F22CCM.area_crop$Et_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F22ccmtrait_info1$Tt_Xmap_Et, ylim = c(0,1),main = "F22 Eup Pheno-> Tetra Pheno", col = "#cfd3db");lines(F22CCM.area_crop$Tt_Xmap_Et[,1], F22CCM.area_crop$Tt_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
# Half 25C and Full 25C
par(mfrow = c(2,2))
plot(rho~LibSize, H25ccmtrait_info1$Et_Xmap_Tt, ylim = c(0,1),main = "H25 Tetra Pheno-> Eup Pheno", col = "#cfd3db");lines(H25CCM.area_crop$Et_Xmap_Tt[,1], H25CCM.area_crop$Et_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, H25ccmtrait_info1$Tt_Xmap_Et, ylim = c(0,1),main = "H25 Eup Pheno-> Tetra Pheno", col = "#cfd3db");lines(H25CCM.area_crop$Tt_Xmap_Et[,1], H25CCM.area_crop$Tt_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
plot(rho~LibSize, F25ccmtrait_info1$Et_Xmap_Tt, ylim = c(0,1),main = "F25 Tetra Pheno-> Eup Pheno", col = "#cfd3db");lines(F25CCM.area_crop$Et_Xmap_Tt[,1], F25CCM.area_crop$Et_Xmap_Tt[,2], col = "#FFC857", lwd = 2)
plot(rho~LibSize, F25ccmtrait_info1$Tt_Xmap_Et, ylim = c(0,1),main = "F25 Eup Pheno-> Tetra Pheno", col = "#cfd3db");lines(F25CCM.area_crop$Tt_Xmap_Et[,1], F25CCM.area_crop$Tt_Xmap_Et[,2], col = "#FFC857", lwd = 2, lty = 2)
5. Appendix II Table S1 Mean CCM Prediction at Largest Library Size by Temperature
Here we create a data frame for all of the highest CCM predictability for each CCM test we ran. We then extracted the mean CCM predictability of each pair at the largest library size (libSizes = 80) to fill in the data frame. “E_Xmap_Tt” in “mapping” column means the CCM test for the casual effects of Tetrahymena phenotypic dynamics on Euplotes ecological dynamics. It’s a bottom-up, written as “B.up”. It tests the effects phenotypic dynamics has on ecological dynamics, thus indicated as “PhenoEco” effects. It is a interaction between Euplotes and Tetrahymena, thus “ET” as species pair. For “pair” column, “T” denotes Tetrahymena, “E” denotes Euplotes, “B” for bacteria,and “Tetra” and “Eup” stands for the effects of one species’ ecological dynamics on its on phenotypic dynamics or vise versa.
# create dataframe to store predictability at largest library size.
area.rho.max <- data.frame(mapping = rep(
c("E_Xmap_Tt", "Tt_Xmap_E", "Et_Xmap_Tt", "Tt_Xmap_Et", "T_Xmap_Tt","Tt_Xmap_T", "Tt_Xmap_B", "B_Xmap_Tt", "E_Xmap_T","T_Xmap_E", "T_Xmap_B", "B_Xmap_T", "E_Xmap_B", "B_Xmap_E","E_Xmap_Et", "Et_Xmap_E", "Et_Xmap_T", "T_Xmap_Et", "Et_Xmap_B","B_Xmap_Et" ),time =4),direc = rep(c("B.up","T.down","B.up","T.down","Single", "Single","B.up","T.down","B.up","T.down","B.up","T.down","B.up","T.down","Single","Single", "B.up","T.down","B.up","T.down"), time = 4), dynamics = rep(c("PhenoEco", "EcoPheno", "PhenoPheno","PhenoPheno","PhenoEco", "EcoPheno","EcoPheno", "PhenoEco", "EcoEco", "EcoEco","EcoEco","EcoEco","EcoEco","EcoEco", "PhenoEco", "EcoPheno","EcoPheno", "PhenoEco","EcoPheno", "PhenoEco"), time = 4),tre = rep(c("Full22C", "Full25C", "Half22C", "Half25C"), each = 20),temp = rep(c("22C", "25C"), each = 20, time = 2),nut =rep(c("Full", "Half"), each = 40),nut2 =rep(c("High", "Low"), each = 40), pair = rep(c("ET","ET","Tetra","BT", "ET", "BT","BE","Eup","ET", "BE"), each = 2, time = 4), rho= rep(NA,80))
# extract predictability at largest library size 80
for (i in 1:8) {
area.rho.max[i,9] <- F22CCM.area_crop[[i]][8,2]
area.rho.max[i+20, 9] <- F25CCM.area_crop[[i]][8,2]
area.rho.max[i+40, 9] <- H22CCM.area_crop[[i]][8,2]
area.rho.max[i+60, 9] <- H25CCM.area_crop[[i]][8,2]
}
for (i in 1:12) {
area.rho.max[i+8,9] <- F22CCM.area[[i]][8,2]
area.rho.max[i+28, 9] <- F25CCM.area[[i]][8,2]
area.rho.max[i+48, 9] <- H22CCM.area[[i]][8,2]
area.rho.max[i+68, 9] <- H25CCM.area[[i]][8,2]
}
head(area.rho.max)## mapping direc dynamics tre temp nut nut2 pair rho
## 1 E_Xmap_Tt B.up PhenoEco Full22C 22C Full High ET 0.9974138
## 2 Tt_Xmap_E T.down EcoPheno Full22C 22C Full High ET 0.7643527
## 3 Et_Xmap_Tt B.up PhenoPheno Full22C 22C Full High ET 0.4844797
## 4 Tt_Xmap_Et T.down PhenoPheno Full22C 22C Full High ET -0.1273202
## 5 T_Xmap_Tt Single PhenoEco Full22C 22C Full High Tetra 0.4272720
## 6 Tt_Xmap_T Single EcoPheno Full22C 22C Full High Tetra 0.2982485Here is the average CCM prediactability accross species pair within each temperature treatment.
# mean rho for temp
rho.mean.temp<- area.rho.max %>% group_by(dynamics, temp) %>% summarise(avg_rho = mean(rho))
rho.mean.temp## # A tibble: 8 × 3
## # Groups: dynamics [4]
## dynamics temp avg_rho
## <chr> <chr> <dbl>
## 1 EcoEco 22C 0.771
## 2 EcoEco 25C 0.743
## 3 EcoPheno 22C 0.551
## 4 EcoPheno 25C 0.657
## 5 PhenoEco 22C 0.707
## 6 PhenoEco 25C 0.786
## 7 PhenoPheno 22C 0.117
## 8 PhenoPheno 25C 0.4546. Appendix III Table S1 All Strong Interaction in the FW
To be conservative and increase the chance of picking up biologically relevant effects, we consider effects whose CCM predictability value was above 0.75 as strong interaction. Here are the 40 strong interactions presented in Figure 5 in the main text the Appendix III Fig S1.
## mapping direc dynamics tre temp nut nut2 pair rho
## 1 E_Xmap_Tt B.up PhenoEco Full22C 22C Full High ET 0.9974138
## 2 Tt_Xmap_E T.down EcoPheno Full22C 22C Full High ET 0.7643527
## 9 E_Xmap_T B.up EcoEco Full22C 22C Full High ET 0.9975146
## 10 T_Xmap_E T.down EcoEco Full22C 22C Full High ET 0.9142912
## 13 E_Xmap_B B.up EcoEco Full22C 22C Full High BE 0.9512350
## 15 E_Xmap_Et Single PhenoEco Full22C 22C Full High Eup 0.9595775
## 19 Et_Xmap_B B.up EcoPheno Full22C 22C Full High BE 0.8551544
## 21 E_Xmap_Tt B.up PhenoEco Full25C 25C Full High ET 0.8514711
## 22 Tt_Xmap_E T.down EcoPheno Full25C 25C Full High ET 0.7576846
## 25 T_Xmap_Tt Single PhenoEco Full25C 25C Full High Tetra 0.9510976
## 26 Tt_Xmap_T Single EcoPheno Full25C 25C Full High Tetra 0.9249959
## 28 B_Xmap_Tt T.down PhenoEco Full25C 25C Full High BT 0.9978469
## 29 E_Xmap_T B.up EcoEco Full25C 25C Full High ET 0.9278876
## 30 T_Xmap_E T.down EcoEco Full25C 25C Full High ET 0.8468363
## 32 B_Xmap_T T.down EcoEco Full25C 25C Full High BT 0.9962212
## 33 E_Xmap_B B.up EcoEco Full25C 25C Full High BE 0.9018809
## 34 B_Xmap_E T.down EcoEco Full25C 25C Full High BE 0.9929382
## 37 Et_Xmap_T B.up EcoPheno Full25C 25C Full High ET 0.8100385
## 40 B_Xmap_Et T.down PhenoEco Full25C 25C Full High BE 0.8913769
## 41 E_Xmap_Tt B.up PhenoEco Half22C 22C Half Low ET 0.9979191
## 42 Tt_Xmap_E T.down EcoPheno Half22C 22C Half Low ET 0.9277864
## 49 E_Xmap_T B.up EcoEco Half22C 22C Half Low ET 0.9885190
## 50 T_Xmap_E T.down EcoEco Half22C 22C Half Low ET 0.8045549
## 52 B_Xmap_T T.down EcoEco Half22C 22C Half Low BT 0.8428261
## 53 E_Xmap_B B.up EcoEco Half22C 22C Half Low BE 0.9686965
## 55 E_Xmap_Et Single PhenoEco Half22C 22C Half Low Eup 0.8837168
## 57 Et_Xmap_T B.up EcoPheno Half22C 22C Half Low ET 0.8909875
## 58 T_Xmap_Et T.down PhenoEco Half22C 22C Half Low ET 0.7527912
## 59 Et_Xmap_B B.up EcoPheno Half22C 22C Half Low BE 0.9587790
## 60 B_Xmap_Et T.down PhenoEco Half22C 22C Half Low BE 0.9255399
## 61 E_Xmap_Tt B.up PhenoEco Half25C 25C Half Low ET 0.9978171
## 62 Tt_Xmap_E T.down EcoPheno Half25C 25C Half Low ET 0.9538906
## 65 T_Xmap_Tt Single PhenoEco Half25C 25C Half Low Tetra 0.9221324
## 69 E_Xmap_T B.up EcoEco Half25C 25C Half Low ET 0.9952086
## 70 T_Xmap_E T.down EcoEco Half25C 25C Half Low ET 0.9350556
## 73 E_Xmap_B B.up EcoEco Half25C 25C Half Low BE 0.7625624
## 75 E_Xmap_Et Single PhenoEco Half25C 25C Half Low Eup 0.9328153
## 77 Et_Xmap_T B.up EcoPheno Half25C 25C Half Low ET 0.7811664
## 78 T_Xmap_Et T.down PhenoEco Half25C 25C Half Low ET 0.7967766
## 79 Et_Xmap_B B.up EcoPheno Half25C 25C Half Low BE 0.8520330